setwd("~/Desktop/grass_exp/")

require(stringr)
require(reshape2)
require(scales)
require(ggplot2)
require(ggthemes)
require(lubridate)
require(ncdf4)
require(tidyr)
require(rstatix)     # needed for adjust_pvalue function
require(ggstar)
library(cowplot)


#source("~/Desktop/tsj_plotting.R")
source("~/Desktop/grass_exp/grass_exp_tools.R")

qe_dir      <- "QE_spectra/"    # directory of raw QE-PRO data files
flm_dir     <- "FLM_spectra/"   # directory of raw FLAME data files
targets_dir <- "targets/"       # directory of all the targeting files

dark_cal_file <- "cal/QEP02110_dark_calibration_raw_values_20200428-1416.csv"   # dark calibration file
rad_cal_file  <- "cal/Rad_Calibration_QEP02110_20210720.csv"                    # radiameteric calibration file
ref_spec_file <- "cal/reference_SIF_spectra.rds"

licor_file <- "licor_data.csv"
pig_file   <- "pigments.csv"


distance <- 18.4  # meters from the optics head to the targeting platform


#met data:-----------------------------------------------------
met_df <- load_NIST_met()
met_df$atm_correction <- compute_atm_factor(met_df$air_temp,
                                            met_df$pressure,
                                            met_df$RH,
                                            distance = distance)

met_df$date <- strftime(met_df$tstamp, format="%Y%m%d")

# data frame organizing the targeting timing:---------------------------------
target_files       <- list.files(targets_dir)
target_array       <- str_split_fixed(target_files, "_",4)
target_df          <- as.data.frame( target_array[,c(2,4)] )
names(target_df)   <- c("date","cycle")
target_df$cycle    <- str_remove(target_df$cycle,".csv")
target_df$qe_file  <- paste0(  qe_dir, "cameron_QEP02110_",  target_df$date , "_cycle_", target_df$cycle, ".nc" )
target_df$flm_file <- paste0( flm_dir, "cameron_FLMS16637_", target_df$date , "_cycle_", target_df$cycle, ".nc" )
target_df$descript <- paste0( target_df$date, "_cycle_", target_df$cycle)

#Load Proximal Measurements:---------------------------------------------------------
sif_df   <- data.frame()
flame_df <- data.frame()
thermal_df <- data.frame()

for(file_idx in 1:length(target_files)){
  cameron_file <- target_df$qe_file[file_idx]
  flame_file   <- target_df$flm_file[file_idx]
  description  <- target_df$descript[file_idx]
  date_str     <- target_df$date[file_idx]
  targets_file <- paste0(targets_dir, target_files[file_idx] )
  targets_df   <- read.csv(targets_file, header = F, stringsAsFactors = F,comment.char = "!")
  names(targets_df) <- c("time","target","surface_temperature")
  targets_df$timestamp <- as.POSIXct( paste0( date_str, targets_df$time), tz="GMT", format = "%Y%m%d%H:%M")
  targets_df$surface_temperature <- (targets_df$surface_temperature - 32)*5/9  #convert F to C


  sif_df <- rbind( sif_df, compute_sif(cameron_file,
                                       targets_df,
                                       dark_cal_file,
                                       rad_cal_file,
                                       ref_spec_file,
                                       met_df) )
  
  flame_df <- rbind( flame_df, compute_flame(flame_file, targets_df))
  tdf <- get_thermal_data( cameron_file, targets_df)
  tdf$cycle = file_idx
  thermal_df <- rbind( thermal_df, tdf )
} 

#thermal check:
thermal_df$date <- date( thermal_df$tstamps )
g <- ggplot( thermal_df, aes(x=cycle,y=spec_temp,group=cycle)) +
  geom_boxplot(outlier.color = NA) +
  facet_wrap(~date,scales="free_x",nrow=1) +
  theme_bw() +
  labs(y="Spectrometer Housing Temperature [°C]")
ggsave("thermal_plot.jpg",g)


ggplot( thermal_df, aes(x=tstamps,y=spec_temp)) +
  geom_point() +facet_wrap(~date,scales="free_x",nrow = 1)

#drop outliers:
sif_df   <- sif_df[ !( (sif_df$description == "20210714_cycle_2") & (sif_df$target == "Drought_4")) ,  ]
sif_df   <- sif_df[ !( (sif_df$description == "20210714_cycle_3") & (sif_df$target == "Drought_4")) ,  ]
flame_df <- flame_df[ !( (flame_df$description == "20210714_cycle_2") & (flame_df$target == "Drought_4")) ,  ]
flame_df <- flame_df[ !( (flame_df$description == "20210714_cycle_3") & (flame_df$target == "Drought_4")) ,  ]


#sif_df$sif_doas <- -10*sif_df$sif_doas

foo <- sif_df[ sif_df$target != "Spect", ]
foo$sif_doas <- foo$sif_doas*2
ggplot(foo, aes(x=sif_fld,sif_doas,color=target_type)) +
  geom_point()

ggplot(foo,aes(x=par)) +
  geom_point(aes(y=sif_doas,color="DOAS")) +
  geom_point(aes(y=sif_fld,color="FLD")) +
  labs(x="PAR",y="SIF",color="")

ggplot(sif_df,aes(x=sif_fld/par,y=sif_doas/par,color=target_type)) +
  geom_point()

# Leaf-Level Licor Data:-----------------------------------------------------------------------
licor_df   <- load_licor_data( licor_file )
licor_df$tstamp <- as.POSIXct( licor_df$tstamp + 4*60*60, origin="1970-01-01", tz="UTC")
licor_df$target_type <- factor( licor_df$target_type,levels=c("Control","ABA","Drought"))

licor_df$long_desc <- "Other"
licor_df$long_desc[ licor_df$target_type == "Drought" & licor_df$date == "20210716"] <- "Drought: Rewatered"
licor_df$long_desc[ licor_df$target_type == "Drought" & licor_df$date == "20210715"] <- "Drought: Day 2"
licor_df$long_desc[ licor_df$target_type == "ABA" & licor_df$date == "20210716"] <- "ABA: Day 3"

# Pigment Data:---------------------------------------------------------------------------------
pig_df <- load_pigments( pig_file )

#aggregate and merge sif and vi data: ------------------------------------------------------
noise_df <- aggregate( cbind(sif_fld,sif_doas,sif_over_par,sif_doas_over_par,par,tstamp,surface_temperature) ~ target+target_type+target_rep+description+date,sif_df,      sd, na.action = NULL )
df_a     <- aggregate( cbind(sif_fld,sif_doas,sif_over_par,sif_doas_over_par,par,tstamp,surface_temperature) ~ target+target_type+target_rep+description+date,sif_df,  median, na.action = NULL )
n_df     <- aggregate( cbind(sif_fld,sif_doas,sif_over_par,sif_doas_over_par,par,tstamp,surface_temperature) ~ target+target_type+target_rep+description+date,sif_df,  length, na.action = NULL )
df_a$sif_fld_sd      <- noise_df$sif_fld
df_a$sif_doas_sd     <- noise_df$sif_doas
df_a$sif_over_par_sd <- noise_df$sif_over_par
df_a$sif_doas_over_par_sd <- noise_df$sif_doas_over_par
df_a$par_sd          <- noise_df$par
df_a$n_sif_samples   <- n_df$sif_fld
df_a$tstamp          <- as.POSIXct(df_a$tstamp, origin="1970-01-01")

noise_df  <- aggregate( cbind(PRI_550,PRI_570,EVI,NDVI,NIRv,tstamp) ~ target+target_type+target_rep+description+date,flame_df,sd, na.action = NULL)
df_b      <- aggregate( cbind(PRI_550,PRI_570,EVI,NDVI,NIRv,tstamp) ~ target+target_type+target_rep+description+date,flame_df,median, na.action = NULL )
n_df      <- aggregate( cbind(PRI_550,PRI_570,EVI,NDVI,NIRv,tstamp) ~ target+target_type+target_rep+description+date,flame_df,length, na.action = NULL )
df_b$PRI_550_sd <- noise_df$PRI_550
df_b$PRI_570_sd <- noise_df$PRI_570
df_b$NDVI_sd    <- noise_df$NDVI
df_b$EVI_sd     <- noise_df$EVI
df_b$NIRv_sd    <- noise_df$NIRv
df_b$n_flame_samples <- n_df$PRI_550
df_b$tstamp   <- as.POSIXct(df_b$tstamp, origin="1970-01-01")
  
merged_df <- merge(df_a,df_b,by=c("target","description","target_type","target_rep","date"))
merged_df$SIF_yield      <- merged_df$sif_over_par/merged_df$NDVI
merged_df$SIF_yield_doas <- merged_df$sif_doas_over_par/merged_df$NDVI
merged_df$SIF_yield_sd  <- merged_df$SIF_yield*(  (merged_df$sif_over_par_sd / merged_df$sif_over_par ) + (merged_df$NDVI_sd / merged_df$NDVI ) )
merged_df$SIF_yield_doas_sd <- merged_df$SIF_yield_doas*( (merged_df$sif_doas_over_par_sd / merged_df$sif_doas_over_par ) + (merged_df$NDVI_sd / merged_df$NDVI ) )

merged_df$desc <- ""
merged_df$desc[ merged_df$date == "20210713"] <- "Day 0"
merged_df$desc[ merged_df$date == "20210714"] <- "Day 1"
merged_df$desc[ merged_df$date == "20210715"] <- "Day 2"
merged_df$desc[ merged_df$date == "20210716"] <- "Day 3"

merged_df$surface_temperature <- replace_na( merged_df$surface_temperature, 0) #Having NAs was messing up the aggregate...just remember day 0 has no surface temps!

small_merged_df <- aggregate( cbind(PRI_550,PRI_570,EVI,NDVI,NIRv,sif_fld,sif_doas,sif_over_par,SIF_yield,SIF_yield_doas,surface_temperature,par)~target_type+target_rep+date,merged_df,"mean")

merged_df$PRI_550_var <- merged_df$PRI_550_sd^2
merged_df$PRI_570_var <- merged_df$PRI_570_sd^2
merged_df$NDVI_var <- merged_df$NDVI_sd^2
merged_df$NIRv_var <- merged_df$NIRv_sd^2
merged_df$EVI_var <- merged_df$EVI_sd^2
merged_df$sif_fld_var <- merged_df$sif_fld_sd^2
merged_df$sif_over_par_var <- merged_df$sif_over_par_sd^2
merged_df$SIF_yield_doas_var <- merged_df$SIF_yield_doas_sd^2
merged_df$SIF_yield_var <- merged_df$SIF_yield_sd^2
merged_df$par_var <- merged_df$par_sd^2

small_merged_noise_df <- aggregate( cbind(PRI_550_var,PRI_570_var,EVI_var,NDVI_var,NIRv_var,sif_fld_var,sif_over_par_var,SIF_yield_var,SIF_yield_doas_var,par_var)~target_type+target_rep+date,merged_df,"sum")
small_merged_n_df <- aggregate( cbind(PRI_550_var)~target_type+target_rep+date,merged_df,"length")
n <- sqrt( small_merged_n_df$PRI_550_var )

small_merged_df$PRI_550_sd <- sqrt( small_merged_noise_df$PRI_550_var ) / n
small_merged_df$PRI_570_sd <- sqrt( small_merged_noise_df$PRI_570_var ) / n
small_merged_df$NDVI_sd    <- sqrt( small_merged_noise_df$NDVI_var ) / n
small_merged_df$NIRv_sd    <- sqrt( small_merged_noise_df$NIRv_var ) / n
small_merged_df$EVI_sd     <- sqrt( small_merged_noise_df$EVI_var ) / n
small_merged_df$sif_fld_sd <- sqrt( small_merged_noise_df$sif_fld_var ) / n
small_merged_df$sif_over_par_sd <- sqrt( small_merged_noise_df$sif_over_par_var ) / n
small_merged_df$SIF_yield_doas_sd <- sqrt( small_merged_noise_df$SIF_yield_doas_var ) / n
small_merged_df$SIF_yield_sd <- sqrt( small_merged_noise_df$SIF_yield_var ) / n
small_merged_df$par_sd <- sqrt( small_merged_noise_df$par_var ) / n

all_df   <- merge(small_merged_df,licor_df,by=c("date","target_type","target_rep"))
all_df   <- merge( all_df, pig_df, by=c("date","target_type","target_rep"))
#all_df <- aggregate( .~date+target_type+target_rep, all_df, "mean", na.action = NULL)
all_df$desc[ all_df$date == "20210713"] <- "Day 0"
all_df$desc[ all_df$date == "20210714"] <- "Day 1"
all_df$desc[ all_df$date == "20210715"] <- "Day 2"
all_df$desc[ all_df$date == "20210716"] <- "Day 3"

all_df$long_desc <- ""
all_df$long_desc[ all_df$target_type == "Drought" & all_df$date == "20210714"] <- "Water Stress"
#all_df$long_desc[ all_df$target_type == "Drought" & all_df$date == "20210715"] <- "Water Stress"
#all_df$long_desc[ all_df$target_type == "Drought" & all_df$date == "20210716"] <- "Water Stress: Rewatered"

all_df$long_desc[ all_df$target_type == "Control" & all_df$date == "20210714"] <- "Control"

all_df$long_desc[ all_df$target_type == "ABA" & all_df$date == "20210714"] <- "ABA"
#all_df$long_desc[ all_df$target_type == "ABA" & all_df$date == "20210715"] <- "ABA"
#all_df$long_desc[ all_df$target_type == "ABA" & all_df$date == "20210716"] <- "ABA"

all_df$alpha_lev <- 1
all_df$alpha_lev[ all_df$long_desc =="" ] <- .7

levs <- c("Day 0 and Controls",
          "ABA: Day 1","ABA: Day 2","ABA: Day 3",
          "Water Stress: Day 1","Water Stress: Day 2","Water Stress: Day 3", "Water Stress: Rewatered")
all_df$long_desc <- factor( all_df$long_desc, levels=levs)

labs(x=bquote(A[net]~~"[ "*mu*mol~m^{-2}*s^{-1}*" ]"), y=bquote(phi[PSII]))

foo <- all_df
levels(foo$target_type) <- c( levels(foo$target_type), "Water Stress")
foo$target_type[ foo$target_type == "Drought"] <- "Water Stress"

ggplot(foo,aes(x=SIF_yield,y=Fs,fill=target_type,starshape=desc)) +
  geom_star(alpha=.8) +
  theme_bw() + 
  scale_starshape_manual(values = c(15,11,13,1)) + 
  scale_fill_manual(values=c("#004488","#BB5566","#DDAA33")) +
  #labs(x=" ", y=" ",title=" ") +
  #theme(legend.position = "none") +
  theme(plot.margin = unit(c(.5, .25, 0, 1), "cm")) 

g1

g1 <- ggplot(foo,aes(x=A,y=PhiPS2,fill=target_type,starshape=desc)) +
  geom_star(alpha=.8) +
  theme_bw() + 
  scale_starshape_manual(values = c(15,11,13,1)) + 
  scale_fill_manual(values=c("#004488","#BB5566","#DDAA33")) +
  labs(x=" ", y=" ",title=" ") +
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(.5, .25, 0, 1), "cm"))
  
g2 <- ggplot(foo,aes(x=A,y=PRI_570,fill=target_type,starshape=desc)) +
  geom_star(alpha=.8) +
  theme_bw() + 
  scale_starshape_manual(values = c(15,11,13,1)) + 
  scale_fill_manual(values=c("#004488","#BB5566","#DDAA33")) +
  labs(x=" ", y=" ",title=" ")+
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(.5, .25, 0, 1), "cm"))

g3 <- ggplot(foo,aes(x=A,y=sif_doas,fill=target_type,starshape=desc)) +
  geom_star(alpha=.8) +
  theme_bw() + 
  scale_starshape_manual(values = c(15,11,13,1)) + 
  scale_fill_manual(values=c("#004488","#BB5566","#DDAA33")) +
  labs(x=" ", y=" ",title=" ") +
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(.5, .25, 0, 1), "cm"))

g4 <- ggplot(foo,aes(x=A,y=1000*SIF_yield,fill=target_type,starshape=desc)) +
  geom_star(alpha=.8) +
  theme_bw() + 
  scale_starshape_manual(values = c(15,11,13,1)) + 
  scale_fill_manual(values=c("#004488","#BB5566","#DDAA33")) +
  labs(x=" ", y=" ",title=" ") +
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(.5, .25, 0, 1), "cm"))

g5 <- ggplot(foo,aes(x=A,y=surface_temperature,fill=target_type,starshape=desc)) +
  geom_star(alpha=.8) +
  theme_bw() + 
  scale_starshape_manual(values = c(15,11,13,1)) + 
  scale_fill_manual(values=c("#004488","#BB5566","#DDAA33")) +
  labs(x=" ", y=" ",title=" ") +
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(.5, .25, 0, 1), "cm")) +
  ylim(18,38)

g <- plot_grid(g1, g2, g3,g4,g5, 
               labels = c('a.', 'b.','c.','d.','e.'), 
               hjust = -1,
               #vjust = 1,
               label_size = 14,
               align = "hv")

g

ggsave("plots/fig3.jpg",g,dpi=300,width=6,height=5,units="in")  




ggplot(foo,aes(x=A,y=PhiPS2,fill=target_type,starshape=desc,xmin=A-A_sd,xmax=A+A_sd,ymin=PhiPS2-PhiPS2_sd,ymax=PhiPS2+PhiPS2_sd)) +
  geom_star(size=3,alpha=.8) +
  geom_errorbarh(alpha=.8,aes(color=target_type)) +
  geom_errorbar(alpha=.8,aes(color=target_type)) +
  theme_bw() +
  scale_shape_manual(values=c(16, 17, 15,18)) + 
  scale_fill_manual(values=c("#004488","#BB5566","#DDAA33","#555555")) +
  scale_color_manual(values=c("#004488","#BB5566","#DDAA33","#555555"),guide="none") +
  scale_starshape_manual(values = c(15,11,13,1)) +
  labs(x=bquote(A[net]~~"[ "*mu*mol~m^{-2}*s^{-1}*" ]"), y=bquote(phi[PSII]),fill="",starshape="")+
  theme(text = element_text(size = 14)) 

ggplot(foo,aes(x=A,y=PRI_570,fill=target_type,starshape=desc)) +
  geom_star(size=3,alpha=.8) +
  theme_bw() +
  scale_shape_manual(values=c(16, 17, 15,18)) + 
  scale_fill_manual(values=c("#004488","#BB5566","#DDAA33")) +
  scale_starshape_manual(values = c(15,11,13,1)) +
  labs(x=bquote(A[net]~~"[ "*mu*mol~m^{-2}*s^{-1}*" ]"), y=bquote(PRI),fill="",starshape="")+
  theme(text = element_text(size = 14))    

ggplot(foo,aes(x=A,y=SIF_yield,fill=target_type,starshape=desc,xmin=A-A_sd,xmax=A+A_sd,ymin=SIF_yield-SIF_yield_sd,ymax=SIF_yield+SIF_yield_sd)) +
  geom_star(size=3,alpha=.8) +
  theme_bw() +
  geom_errorbarh(alpha=.7,aes(color=target_type)) +
  geom_errorbar(alpha=.7,aes(color=target_type)) +
  scale_shape_manual(values=c(16, 17, 15,18)) + 
  scale_fill_manual(values=c("#004488","#BB5566","#DDAA33")) +
  scale_color_manual(values=c("#004488","#BB5566","#DDAA33","#555555"),guide="none") +
  scale_starshape_manual(values = c(15,11,13,1)) +
  labs(x=bquote(A[net]~~"[ "*mu*mol~m^{-2}*s^{-1}*" ]"), y=bquote(SIF[yield]),fill="",starshape="") +
  theme(text = element_text(size = 14))    

ggplot(foo,aes(x=FvFm,y=PRI_570,fill=target_type,starshape=desc)) +
  geom_star(size=3,alpha=.8) +
  theme_bw() +
  #geom_errorbarh(alpha=.7,aes(color=target_type)) +
  #geom_errorbar(alpha=.7,aes(color=target_type)) +
  scale_shape_manual(values=c(16, 17, 15,18)) + 
  scale_fill_manual(values=c("#004488","#BB5566","#DDAA33")) +
  scale_color_manual(values=c("#004488","#BB5566","#DDAA33","#555555"),guide="none") +
  scale_starshape_manual(values = c(15,11,13,1)) +
  labs(x=bquote(FvFm), y=bquote(PRI_570),fill="",starshape="") +
  theme(text = element_text(size = 14))  

ggplot(foo,aes(x=SIF_yield,y=PRI_570,fill=target_type,starshape=desc)) +
  geom_star(size=3,alpha=.8) +
  theme_bw() +
  scale_shape_manual(values=c(16, 17, 15,18)) + 
  scale_fill_manual(values=c("#004488","#BB5566","#DDAA33")) +
  scale_starshape_manual(values = c(15,11,13,1)) +
  labs(x=bquote(SIF[yield]), y=bquote(PRI),fill="",starshape="") +
  theme(text = element_text(size = 14))   
#Figure 3:----------------------
#a:
gdf <- data.frame()
for(date in unique(licor_df$date)){
  ldf <- licor_df
  ldf$highlight <- date
  ldf$long_desc <- ""
  ldf$long_desc[ ldf$date == date ] <- as.character( ldf$target_type[ ldf$date == date] )
  gdf <- rbind( gdf, ldf)
}

gdf$highlight[ gdf$highlight == "20210713" ] <- "Day 0"
gdf$highlight[ gdf$highlight == "20210714" ] <- "Day 1"
gdf$highlight[ gdf$highlight == "20210715" ] <- "Day 2"
gdf$highlight[ gdf$highlight == "20210716" ] <- "Day 3"

gdf$alph <- 1
gdf$alph[ gdf$long_desc == ""] <- .3

gdf$fit_data <- gdf$PhiPS2
gdf$fit_data[ gdf$long_desc == ""] <- NA 

gdf$long_desc <- factor( gdf$long_desc, levels= c("","ABA","Control","Drought","Water Stress"))
gdf$long_desc[ gdf$long_desc == "Drought"] <- "Water Stress"

g <- ggplot(gdf,aes(x=A,y=PhiPS2,xmin=A-A_sd,xmax=A+A_sd,ymin=PhiPS2-PhiPS2_sd,ymax=PhiPS2+PhiPS2_sd,color=long_desc)) +
  geom_point(stroke=0,size=2.5,alpha=gdf$alph) +
  geom_errorbarh(alpha=gdf$alph) +
  geom_errorbar(alpha=gdf$alph) +
  theme_bw() +
  scale_color_manual(values=c("#000000","#B51D14","#4053D3","#DDB310")) +
  #scale_shape_manual(values=c(3, 16, 17)) +
  labs(x=bquote(A[net]~~"[ "*mu*mol~m^{-2}*s^{-1}*" ]"), y=bquote(SIF[yield])) +
  geom_smooth(aes(group=highlight), se=F,method="lm",fullrange=T,linetype=2,color="black",size=.5) +
  #geom_smooth(aes(group=highlight,y=fit_data), se=F,method="lm",fullrange=T) +
  facet_wrap(~highlight, nrow=1) 
  #theme(legend.position = "none")
ggsave("plots/fig3_part1.jpg",g,dpi=300,width=8,height=3,units="in")


ggplot( all_df,aes(x=PRI_570,y=sif_fld,color=target_type)) +
  geom_point() 
  #facet_wrap(~date)
ggplot( all_df,aes(x=A,y=SIF_yield,color=target_type,group=paste0(target_type,target_rep))) +
  geom_point()

foo <- aggregate(cbind(A,PRI_570,FvFm,PhiPS2,SIF_yield)~date+target_type,all_df,"mean")

ggplot( foo2,aes(x=A,y=SIF_yield,color=target_type) ) +
  geom_point() +
  geom_path(arrow=arrow())+
  geom_point(data=foo,alpha=.7)

ggplot( foo,aes(x=A,y=PRI_570,color=target_type) ) +
  geom_point() +
  geom_path(arrow=arrow())

ggplot( foo2,aes(x=SIF_yield,y=PRI_570,color=target_type) ) +
  geom_point() +
  geom_path(arrow=arrow())+
  geom_point(data=foo,alpha=.5)

ggplot( foo2,aes(x=A,y=PhiPS2,color=target_type) ) +
  geom_point() +
  geom_path(arrow=arrow())+
  geom_point(data=foo,alpha=.7)


ggplot( foo2,aes(x=A,y=FvFm,color=target_type) ) +
  geom_point() +
  geom_path(arrow=arrow()) +
  geom_point(data=foo,alpha=.7)

foo <- all_df[ all_df$date %in% c(20210713,20210714,20210715), ]
foo2 <- aggregate(cbind(A,PRI_570,FvFm,PhiPS2,SIF_yield)~date+target_type,foo,"mean")
ggplot( foo2,aes(x=A,y=FvFm,color=target_type) ) +
  geom_path(arrow=arrow()) +
  geom_point(data=foo,alpha=.7)
names(foo)

#b:
gdf <- data.frame()
for(date in unique(df_a$date)){
  ldf <- df_a
  ldf$highlight <- date
  ldf$long_desc <- ""
  ldf$long_desc[ ldf$date == date ] <- as.character( ldf$target_type[ ldf$date == date] )
  gdf <- rbind( gdf, ldf)
}

gdf$highlight[ gdf$highlight == "20210713" ] <- "Day 0"
gdf$highlight[ gdf$highlight == "20210714" ] <- "Day 1"
gdf$highlight[ gdf$highlight == "20210715" ] <- "Day 2"
gdf$highlight[ gdf$highlight == "20210716" ] <- "Day 3"

gdf$alph <- .8
gdf$alph[ gdf$long_desc == ""] <- .2

gdf$fit_data <- gdf$sif_fld
gdf$fit_data[ gdf$long_desc == ""] <- NA 

gdf$long_desc <- factor( gdf$long_desc, levels= c("","ABA","Control","Drought","Water Stress"))
gdf$long_desc[ gdf$long_desc == "Drought"] <- "Water Stress"

ggplot(gdf,aes(x=par,y=sif_fld,color=long_desc,xmin=par-par_sd,xmax=par+par_sd,ymin=sif_fld-sif_fld_sd,ymax=sif_fld+sif_fld_sd)) +
  geom_point(stroke=0,size=2,alpha=gdf$alph) +
  #geom_path(aes(group=paste0(target_rep,target_type)))+
  geom_errorbarh(alpha=gdf$alph) +
  geom_errorbar(alpha=gdf$alph) +
  scale_color_manual(values=c("#000000","#B51D14","#4053D3","#DDB310")) +  theme_bw()+
  labs(x=bquote(PAR~~"[ "*mu*mol~m^{-2}*s^{-1}*" ]"), y=bquote(SIF~~"[ "*mW~m^{-2}*s^{-1}*sr^{-1}*nm^{-1}*" ]")) +
  geom_smooth(aes(group=highlight), se=F,method="lm",fullrange=T,linetype=2,color="black",size=.5) +
  #geom_smooth(aes(group=highlight,y=fit_data), se=F,method="lm",fullrange=T) +
  facet_wrap(~highlight, nrow=1) +
  theme(legend.position = "none")

g

ggsave("plots/fig3_part2_allpoints.jpg",g,dpi=300,width=8,height=3,units="in")


#percent difference bootstrap:-------------------------
df <- melt( all_df, id=c("date","target_type","target_rep"))
#df <- melt( sif_df, id=c("description","target_type","target_rep"))
t_df <- data.frame()
dates <- unique(df$date)
descriptions <- unique(df$description)
target_reps <- unique(df$target_rep)
vars <- c("sif_over_par", "NIRv","PRI_570", "A","gsw","NDVI","sif_fld","sif_doas","sif_over_par","SIF_yield","SIF_yield_doas","FvFm","surface_temperature","Fmp","Fs") #,"C_pool","Cla_frac","X_pool","ZA_per_X")
#names(sds) <- vars
n_repeat <- 1000
n_replicates <- 5
set.seed(5)
use_error = 0
for( date in dates){
  #for( desc in descriptions){
  print(date)
  for( var in vars){
    print(var)
    #sample_list  <- sample(target_reps,n_repeat, replace=TRUE)
    #control_err <- rnorm( n_replicates*n_repeat,mean=0,sd=sds[var]) 
    #sample_err  <- rnorm( n_replicates*n_repeat,mean=0,sd=sds[var])
    for( target_type in c("ABA","Drought")){
      
      c_df <- df[ (df$date == date) & 
                    (df$target_type == "Control") & 
                    (df$variable == var), ]
      
      c_err_df <- df[ (df$date == date) & 
                        (df$target_type == "Control") & 
                        (df$variable == paste0(var,"_sd")), ]
      
      a_df <- df[   (df$date == date) & 
                      (df$target_type == target_type) & 
                      (df$variable == var), ]
      
      mean_diff <-  (mean(as.numeric(a_df$value)) - mean(as.numeric(c_df$value)) ) / mean(as.numeric(c_df$value))
      
      if( is.na(c_df$value) ){
        print("No data found!")
      }else{
        c_samples <- sample(nrow(c_df), n_replicates*n_repeat, replace=TRUE)
        c_df <- c_df[c_samples, ]
        c_df$repeat_number <- rep( 1:n_repeat, n_replicates)
        c_err_df <- c_err_df[c_samples, ]
        
        a_df <- df[   (df$date == date) & 
                        (df$target_type == target_type) & 
                        (df$variable == var), ]
        
        a_err_df <- df[ (df$date == date) & 
                          (df$target_type == target_type) & 
                          (df$variable == paste0(var,"_sd")), ]
        
        a_samples <- sample(nrow(a_df), n_replicates*n_repeat, replace=TRUE)
        a_df <- a_df[a_samples, ]
        a_df$repeat_number <- rep( 1:n_repeat, n_replicates)
        a_err_df <- a_err_df[a_samples, ]
        a_err <-c()
        c_err <-c()
        if(use_error){
          for(ii in 1:length(c_err_df$value)){
            c_err    <- c(c_err,rnorm( 1, mean=0,sd= as.numeric( c_err_df$value[ii]) ) )
            a_err    <- c(a_err,rnorm( 1, mean=0, sd= as.numeric(a_err_df$value[ii])) )
          }
        }else{
          a_err <- 0
          c_err <- 0
        }
        
        a_df$a_err <- a_err
        c_df$c_err <- c_err
        a_df$value <- as.numeric( a_df$value)
        c_df$value <- as.numeric( c_df$value)
        a_df <- aggregate( value~date+repeat_number, a_df,"mean")
        c_df <- aggregate( value~date+repeat_number, c_df,"mean")
        
        vals_frac  <- ( (as.numeric(a_df$value) + a_err ) / (as.numeric(c_df$value) + c_err)) - 1
        vals_unit  <- ( (as.numeric(a_df$value) + a_err ) - (as.numeric(c_df$value) + c_err))
        t_df <- rbind( t_df , data.frame( date = date,
                                          target_type = target_type,
                                          var = var, 
                                          val_frac = vals_frac,
                                          val_unit = vals_unit))
      }
    }
  }
}


t_df$desc <- ""
t_df$desc[ t_df$date == "20210713"] <- "B.T."
t_df$desc[ t_df$date == "20210714"] <- "day 1"
t_df$desc[ t_df$date == "20210715"] <- "day 2"
t_df$desc[ t_df$date == "20210716"] <- "day 3"

#Remove the drought day 3 for plotting...if you want.
#t_df$val_frac[ (t_df$long_desc == "C4") ] <- NA
#t_df$vals_unit[ (t_df$long_desc == "C4") ] <- NA


#FIGURE 2:----------------------------------------------------------------------------------------------
mdf    <- aggregate(cbind(val_frac,val_unit)~target_type+var+date,t_df,"median")
lo_df  <- aggregate(cbind(val_frac,val_unit)~target_type+var+date,t_df,FUN=quantile, probs=.025)
hi_df  <- aggregate(cbind(val_frac,val_unit)~target_type+var+date,t_df,FUN=quantile, probs=.975)
mdf$lo_frac <- lo_df$val_frac
mdf$hi_frac <- hi_df$val_frac
mdf$lo_unit <- lo_df$val_unit
mdf$hi_unit <- hi_df$val_unit

vars <- c("gsw","A","FvFm","Fmp","Fs","PhiPS2","PRI_570","NDVI","sif_fld","sif_doas","sif_over_par","SIF_yield","SIF_yield_doas", "NIRv","surface_temperature")
formatted_labs <- c("gsw","A[net]","F[v]*'′'*'/'*F[m]*'′'","Fmp","Fs","PHIPSII","PRI","NDVI","SIF","SIFdoas","SIFoverPAR","SIF[yield]","SIF[yield]DOAS","NIRv","T[surf]")
formatted_labs <- paste0("italic(",formatted_labs,")")


# T-TEST FOR SIGNIFICANCE:-----------------------------------------------------------------------------
out_df <- data.frame()
for(var in vars){
  for(t in c("ABA","Drought")){
    for(date in dates){
    sub_df <- df[ (df$variable==var) & (df$target_type %in% c(t,"Control")) & (df$date == date), ]
    sub_df$value <- as.numeric( sub_df$value )
    
    one.way_adj <- sub_df %>%
      anova_test(dv = value, wid = target_rep, within = target_type) %>%
      get_anova_table() %>%
      adjust_pvalue(method = "bonferroni")
    
    ttest <- pairwise.t.test(sub_df$value,
                            sub_df$target_type,
                            paired = F,
                            p.adjust.method = "bonferroni")
    
    out_df <- rbind(out_df, data.frame(var=var,target_type=t, date=date, p=as.numeric(ttest$p.value) ) )
    }
  }
}
out_df$stars <- ""
out_df$stars[ out_df$p < .05 ] <- "*"
out_df$stars[ out_df$p < .01 ] <- "**"
out_df$stars[ out_df$p < .001 ] <- "***"


sub_df <- mdf[ mdf$var %in% vars, ]
sub_df <- merge(sub_df,out_df,by=c("date","target_type","var"))
sub_df$factor <- factor( sub_df$var, levels=vars,labels=formatted_labs)

sub_df$day <- 0
sub_df$day[ sub_df$date == "20210714"] <- 1
sub_df$day[ sub_df$date == "20210715"] <- 2
sub_df$day[ sub_df$date == "20210716"] <- 3
levels( sub_df$target_type ) <- c("ABA","Drought","Water Stress")
sub_df$target_type[ sub_df$target_type =="Drought"] <- "Water Stress"


sdf <- sub_df[ sub_df$var == "gsw", ]
g1 <- ggplot(sdf,aes(x=day,y=val_frac*100,ymin=lo_frac*100,ymax=hi_frac*100,label="*")) +
  geom_abline( slope=0,intercept = 0, size=.3) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(position = position_dodge(width = 0.4),width=.3,size=.3) +
  geom_text(aes(y=Inf,vjust=1.5, label=stars),size=5,family="sans") +
  theme_bw() +
  facet_wrap(~target_type,ncol=2) +
  labs(x="",y="",title = "") +
  theme(legend.position = "none") +
  geom_vline(data=data.frame(xint=0.6, target_type="ABA"), aes(xintercept=xint) ,size=.7,linetype=2,color="red") +
  geom_vline(data=data.frame(xint=2.6, target_type="Water Stress"), aes(xintercept=xint) ,size=.7,linetype=2,color="blue") +
  theme( strip.text = element_text(size = 12),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 17),
         plot.title = element_text(hjust = -.1, size=20) ) +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
#ggsave("plots/fig2a.jpg",g,dpi=300,width=4,height=3,units="in")

g1
sdf <- sub_df[ sub_df$var == "A", ]
g2 <- ggplot(sdf,aes(x=day,y=val_frac*100,ymin=lo_frac*100,ymax=hi_frac*100,label="*")) +
  geom_abline( slope=0,intercept = 0, size=.3) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(position = position_dodge(width = 0.4),width=.3,size=.3) +
  geom_text(aes(y=Inf,vjust=1.5, label=stars),size=5,family="sans") +
  theme_bw() +
  #geom_text(hjust=0, vjust=-5,position = position_dodge(width = 0.8))+
  facet_wrap(~target_type,ncol=2) +
  labs(x="",y="", title = "") +
  theme(legend.position = "right") +
  geom_vline(data=data.frame(xint=0.6, target_type="ABA"), aes(xintercept=xint) ,size=.7,linetype=2,color="red") +
  geom_vline(data=data.frame(xint=2.6, target_type="Water Stress"), aes(xintercept=xint) ,size=.7,linetype=2,color="blue") +
  theme( strip.text = element_text(size = 12),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 17),
         plot.title = element_text(hjust = -.1, size=20) )+
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
#ggsave("plots/fig2b.jpg",g,dpi=300,width=4,height=3,units="in")

sdf <- sub_df[ sub_df$var == "FvFm", ]
g3 <- ggplot(sdf,aes(x=day,y=val_frac*100,ymin=lo_frac*100,ymax=hi_frac*100,label="*")) +
  geom_abline( slope=0,intercept = 0, size=.3) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(position = position_dodge(width = 0.4),width=.3,size=.3) +
  geom_text(aes(y=Inf,vjust=1.5, label=stars),size=5,family="sans") +
  theme_bw() +
  #geom_text(hjust=0, vjust=-5,position = position_dodge(width = 0.8))+
  facet_wrap(~target_type,ncol=2) +
  labs(x="",y="", title = "") +
  theme(legend.position = "right") +
  geom_vline(data=data.frame(xint=0.6, target_type="ABA"), aes(xintercept=xint) ,size=.7,linetype=2,color="red") +
  geom_vline(data=data.frame(xint=2.6, target_type="Water Stress"), aes(xintercept=xint) ,size=.7,linetype=2,color="blue") +
  theme( strip.text = element_text(size = 12),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 17))+
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
#ggsave("plots/fig2c.jpg",g,dpi=300,width=4,height=3,units="in")

sdf <- sub_df[ sub_df$var == "NDVI", ]
g4 <- ggplot(sdf,aes(x=day,y=val_frac*100,ymin=lo_frac*100,ymax=hi_frac*100,label="*")) +
  geom_abline( slope=0,intercept = 0, size=.3) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(position = position_dodge(width = 0.4),width=.3,size=.3) +
  geom_text(aes(y=Inf,vjust=1.5, label=stars),size=5,family="sans") +
  theme_bw() +
  #geom_text(hjust=0, vjust=-5,position = position_dodge(width = 0.8))+
  facet_wrap(~target_type,ncol=2) +
  labs(x="",y="", title = "" ) +
  theme(legend.position = "right") +
  geom_vline(data=data.frame(xint=0.6, target_type="ABA"), aes(xintercept=xint) ,size=.7,linetype=2,color="red") +
  geom_vline(data=data.frame(xint=2.6, target_type="Water Stress"), aes(xintercept=xint) ,size=.7,linetype=2,color="blue") +
  theme( strip.text = element_text(size = 12),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 17) )+
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
#ggsave("plots/fig2d.jpg",g,dpi=300,width=4,height=3,units="in")

sdf <- sub_df[ sub_df$var == "NIRv", ]
g5 <- ggplot(sdf,aes(x=day,y=val_frac*100,ymin=lo_frac*100,ymax=hi_frac*100,label="*")) +
  geom_abline( slope=0,intercept = 0, size=.3) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(position = position_dodge(width = 0.4),width=.3,size=.3) +
  geom_text(aes(y=Inf,vjust=1.5, label=stars),size=5,family="sans") +
  theme_bw() +
  #geom_text(hjust=0, vjust=-5,position = position_dodge(width = 0.8))+
  facet_wrap(~target_type,ncol=2) +
  labs(x="",y="", title = "" ) +
  theme(legend.position = "right") +
  geom_vline(data=data.frame(xint=0.6, target_type="ABA"), aes(xintercept=xint) ,size=.7,linetype=2,color="red") +
  geom_vline(data=data.frame(xint=2.6, target_type="Water Stress"), aes(xintercept=xint) ,size=.7,linetype=2,color="blue") +
  theme( strip.text = element_text(size = 12),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 17) ) +
  ylim(-40,25)+
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
#ggsave("plots/fig2e.jpg",g,dpi=300,width=4,height=3,units="in")

sdf <- sub_df[ sub_df$var == "PRI_570", ]
g6 <- ggplot(sdf,aes(x=day,y=val_frac*100,ymin=lo_frac*100,ymax=hi_frac*100,label="*")) +
  geom_abline( slope=0,intercept = 0, size=.3) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(position = position_dodge(width = 0.4),width=.3,size=.3) +
  geom_text(aes(y=Inf,vjust=1.5, label=stars),size=5,family="sans") +
  theme_bw() +
  #geom_text(hjust=0, vjust=-5,position = position_dodge(width = 0.8))+
  facet_wrap(~target_type,ncol=2) +
  labs(x="",y="", title = "" ) +
  theme(legend.position = "right") +
  geom_vline(data=data.frame(xint=0.6, target_type="ABA"), aes(xintercept=xint) ,size=.7,linetype=2,color="red") +
  geom_vline(data=data.frame(xint=2.6, target_type="Water Stress"), aes(xintercept=xint) ,size=.7,linetype=2,color="blue") +
  theme( strip.text = element_text(size = 12),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 17) )+
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
#ggsave("plots/fig2f.jpg",g,dpi=300,width=4,height=3,units="in")

sdf <- sub_df[ sub_df$var == "sif_fld", ]
g7 <- ggplot(sdf,aes(x=day,y=val_frac*100,ymin=lo_frac*100,ymax=hi_frac*100,label="*")) +
  geom_abline( slope=0,intercept = 0, size=.3) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(position = position_dodge(width = 0.4),width=.3,size=.3) +
  geom_text(aes(y=Inf,vjust=1.5, label=stars),size=5,family="sans") +
  theme_bw() +
  #geom_text(hjust=0, vjust=-5,position = position_dodge(width = 0.8))+
  facet_wrap(~target_type,ncol=2) +
  labs(x="",y="", title = "") +
  theme(legend.position = "right") +
  geom_vline(data=data.frame(xint=0.6, target_type="ABA"), aes(xintercept=xint) ,size=.7,linetype=2,color="red") +
  geom_vline(data=data.frame(xint=2.6, target_type="Water Stress"), aes(xintercept=xint) ,size=.7,linetype=2,color="blue") +
  theme( strip.text = element_text(size = 12),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 17) )+
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
#ggsave("plots/fig2g.jpg",g,dpi=300,width=4,height=3,units="in")

sdf <- sub_df[ sub_df$var == "SIF_yield", ]
g8 <- ggplot(sdf,aes(x=day,y=val_frac*100,ymin=lo_frac*100,ymax=hi_frac*100,label="*")) +
  geom_abline( slope=0,intercept = 0, size=.3) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(position = position_dodge(width = 0.4),width=.3,size=.3) +
  geom_text( aes(y=Inf,vjust=1.5, label=stars),size=5,family="sans") +
  theme_bw() +
  #geom_text(hjust=0, vjust=-5,position = position_dodge(width = 0.8))+
  facet_wrap(~target_type,ncol=2) +
  labs(x="",y="", title = "" ) +
  theme(legend.position = "none") +
  geom_vline(data=data.frame(xint=0.6, target_type="ABA"), aes(xintercept=xint) ,size=.7,linetype=2,color="red") +
  geom_vline(data=data.frame(xint=2.6, target_type="Water Stress"), aes(xintercept=xint) ,size=.7,linetype=2,color="blue") +
  theme( strip.text = element_text(size = 12),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 17))+
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))

g <- plot_grid(g1, g2, g3,g4,g5,g6,g7,g8, 
          labels = c('a.', 'b.','c.','d.','e.','f.','g.','h.'), 
          hjust = -2,
          vjust = 1.5,
          label_size = 14,
          align = "hv")

ggsave("plots/fig2.jpg",g,dpi=300,width=10,height=8,units="in")


sdf <- sub_df[ sub_df$var == "Fs", ]
ggplot(sdf,aes(x=day,y=val_frac*100,ymin=lo_frac*100,ymax=hi_frac*100,label="*")) +
  geom_abline( slope=0,intercept = 0, size=.3) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(position = position_dodge(width = 0.4),width=.3,size=.3) +
  geom_text(aes(y=Inf,vjust=1.5, label=stars),size=5,family="sans") +
  theme_bw() +
  facet_wrap(~target_type,ncol=2) +
  labs(x="",y="",title = "") +
  theme(legend.position = "none") +
  geom_vline(data=data.frame(xint=0.6, target_type="ABA"), aes(xintercept=xint) ,size=.7,linetype=2,color="red") +
  geom_vline(data=data.frame(xint=2.6, target_type="Water Stress"), aes(xintercept=xint) ,size=.7,linetype=2,color="blue") +
  theme( strip.text = element_text(size = 12),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 17),
         plot.title = element_text(hjust = -.1, size=20) ) +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))

#sdf <- sub_df[ sub_df$var %in% c("sif_fld","sif_doas"), ]
sdf <- sub_df[ sub_df$var %in% c("SIF_yield","SIF_yield_doas"), ]
g8 <- ggplot(sdf,aes(x=day,y=val_frac*100,ymin=lo_frac*100,ymax=hi_frac*100,label="*",color=var)) +
  geom_abline( slope=0,intercept = 0, size=.3) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(position = position_dodge(width = 0.4),width=.3,size=.3) +
  geom_text( aes(y=Inf,hjust=1.25, label=stars),position = position_dodge(width = 0.4), size=5,family="sans",angle=90, show.legend = FALSE) +
  theme_bw() +
  #geom_text(hjust=0, vjust=-5,position = position_dodge(width = 0.8))+
  facet_wrap(~target_type,ncol=2) +
  labs(x="Day of Experiment",y="% Difference From Control", color = "Method" ) +
  #theme(legend.position = "none") +
  geom_vline(data=data.frame(xint=0.6, target_type="ABA"), aes(xintercept=xint) ,size=.7,linetype=2,color="red") +
  geom_vline(data=data.frame(xint=2.6, target_type="Water Stress"), aes(xintercept=xint) ,size=.7,linetype=2,color="blue") +
  theme( strip.text = element_text(size = 12),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 14)) +
  scale_color_manual(values = c("#000000","#22BBEE"))
g8

sdf <- sub_df[ sub_df$var == "surface_temperature", ]
levels( sdf$target_type ) <- c("ABA","Drought","Water Stress")
sdf$target_type[ sdf$target_type =="Drought"] <- "Water Stress"
g <- ggplot(sdf,aes(x=day,y=val_unit,ymin=lo_unit,ymax=hi_unit,label="*")) +
  geom_abline( slope=0,intercept = 0, size=.3) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(position = position_dodge(width = 0.4),width=.3,size=.3) +
  geom_text(aes(y=Inf,vjust=1.5, label=stars),size=5,family="sans") +
  theme_bw() +
  #geom_text(hjust=0, vjust=-5,position = position_dodge(width = 0.8))+
  facet_wrap(~target_type,ncol=2) +
  labs(x="",y="", title = bquote( italic( "Surface Temperature")~ "[°C]"  ) ) +
  theme(legend.position = "right") +
  geom_vline(data=data.frame(xint=0.6, target_type="ABA"), aes(xintercept=xint) ,size=.7,linetype=2,color="red") +
  geom_vline(data=data.frame(xint=2.6, target_type="Water Stress"), aes(xintercept=xint) ,size=.7,linetype=2,color="blue") +
  theme( strip.text = element_text(size = 12),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 17),
         plot.title = element_text(hjust = 0.5, size = 18, family = "serif")) +
  expand_limits(y = 10) +
  xlim(0,3.25)
g

lm_eqn <- function(df){
  m <- lm(sif_doas ~ sif_fld, df);
  eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
                   list(a = format(unname(coef(m)[1]), digits = 2),
                        b = format(unname(coef(m)[2]), digits = 2),
                        r2 = format(summary(m)$r.squared, digits = 3)))
  as.character(as.expression(eq));
}

ggplot(sif_df,aes(x=sif_fld,y=sif_doas)) +
  geom_point(alpha=.3) +
  theme_bw() +
  geom_abline(intercept = 0,slope=1,color="red") +
  geom_smooth(method = "lm", se=FALSE, color="blue", formula = y ~ x) 

foo <- lm( sif_doas ~ sif_fld, sif_df )
sum
##---------------------------------
sdf <- sub_df[ sub_df$var %in% c("A","SIF_yield","PRI_570"), ]
ggplot(sdf,aes(x=day,y=val_frac,fill=var)) +
  geom_bar(position = "stack",stat="identity") +
  facet_wrap(~target_type)


ggplot(all_df,aes(x=PhiPS2,y=FvFm,color=SIF_yield)) +
  geom_point() +
  scale_color_viridis_c(limits = c(-3,3)) +
  xlim( 0,.2) +
  ylim(.15,3)

ggplot(licor_df,aes(x=PhiPS2,FvFm,color=target_type)) +
  geom_point() +
  facet_wrap(~date)

#ORIGINAL VERSIONS OF PLOT 3---------------------------------------------------------------

daily_pal <- c('#EEEEEE','#E7D4E8','#C2A5CF', '#9970AB',
               '#FEE391', '#FEC44F', '#81C4E7')



all_df$long_desc <- "A"
all_df$long_desc[ all_df$date == "20210713"] <- "A"
all_df$long_desc[ all_df$date == "20210714"&(all_df$target_type == "ABA")] <- "B1"
all_df$long_desc[ all_df$date == "20210715"&(all_df$target_type == "ABA")] <- "B2"
all_df$long_desc[ all_df$date == "20210716"&(all_df$target_type == "ABA")] <- "B3"
all_df$long_desc[ (all_df$date == "20210716")&(all_df$target_type == "ABA")] <- "B4"
all_df$long_desc[ all_df$date == "20210714"&(all_df$target_type == "Drought")] <- "C1"
all_df$long_desc[ all_df$date == "20210715"&(all_df$target_type == "Drought")] <- "C2"
all_df$long_desc[ all_df$date == "20210716"&(all_df$target_type == "Drought")] <- "C3"
all_df$long_desc[ (all_df$date == "20210716")&(all_df$target_type == "Drought")] <- "C4"

merged_df$long_desc <- "A"
merged_df$long_desc[ merged_df$date == "20210713"] <- "A"
merged_df$long_desc[ merged_df$date == "20210714"&(merged_df$target_type == "ABA")] <- "B1"
merged_df$long_desc[ merged_df$date == "20210715"&(merged_df$target_type == "ABA")] <- "B2"
merged_df$long_desc[ merged_df$date == "20210716"&(merged_df$target_type == "ABA")] <- "B3"
merged_df$long_desc[ (merged_df$date == "20210716")&(merged_df$target_type == "ABA")] <- "B4"
merged_df$long_desc[ merged_df$date == "20210714"&(merged_df$target_type == "Drought")] <- "C1"
merged_df$long_desc[ merged_df$date == "20210715"&(merged_df$target_type == "Drought")] <- "C2"
merged_df$long_desc[ merged_df$date == "20210716"&(merged_df$target_type == "Drought")] <- "C3"
merged_df$long_desc[ (merged_df$date == "20210716")&(merged_df$target_type == "Drought")] <- "C4"

ggplot(all_df,aes(x=A,y=SIF_yield,color=long_desc)) +
  geom_point() + 
  #scale_color_manual(values=daily_pal) +
  #theme_tsj() +
  geom_errorbarh(aes(xmin=A-A_sd,xmax=A+A_sd)) +
  geom_errorbar(aes(ymin=(SIF_yield-SIF_yield_sd),ymax=(SIF_yield+SIF_yield_sd))) +
  scale_color_manual(values=c('#999999','#E7D4E8','#C2A5CF', '#9970AB',
                              '#FEE391', '#FEC44F', '#81C4E7')) +
  #scale_color_manual(values = c("#BB5566","#DDAA33","#004488","#888888")) +
  #scale_color_manual(values = c('#999999', '#FB9A29', '#EC7014', '#993404') ) +
  labs(y=bquote(SIF[yield]~~~"["*A.U.*"]"), x=bquote(A[net]~~~'['*mu*mol~m^-2*s^-1*']')) +
  theme(legend.position = "none")

fit <- lm(SIF_yield~A,all_df,weights = 1/(SIF_yield_sd^2 + A_sd^2) )
summary(fit)

fit <- lm(SIF_yield2~A,all_df,weights = 1/(SIF_yield2_sd^2 + A_sd^2))
summary(fit)

fit <- lm(sif_over_par~A,all_df,weights = 1/(sif_over_par_sd^2 + A_sd^2))
summary(fit)

ggplot(all_df,aes(x=A,y=sif_fld,color=long_desc)) +
  geom_point() + 
  scale_color_manual(values=c('#999999','#E7D4E8','#C2A5CF', '#9970AB',
                              '#FEE391', '#FEC44F', '#81C4E7')) +
  geom_errorbarh(aes(xmin=A-A_sd,xmax=A+A_sd)) +
  #geom_errorbar(aes(ymin=(sif_over_par-sif_over_par_sd)/4600,ymax=(sif_over_par+sif_over_par_sd)/4600)) +
  #scale_color_manual(values = c("#BB5566","#DDAA33","#004488","#888888")) +
  #scale_color_manual(values = c('#999999', '#FB9A29', '#EC7014', '#993404') ) +
  labs(y=bquote("SIF"~~~"["*sr^-1*nm^-1~"]"), x=bquote(A[net]~~~'['*mu*mol~m^-2*s^-1*']'))

ggplot(all_df,aes(x=C_pool,y=X_pool,color=long_desc)) +
  geom_point(size=2) + 
  scale_color_manual(values=c('#999999','#E7D4E8','#C2A5CF', '#9970AB',
                              '#FEE391', '#FEC44F', '#81C4E7')) +
  theme_tsj() +
  #geom_errorbarh(aes(xmin=A-A_sd,xmax=A+A_sd)) +
  #geom_errorbar(aes(ymin=(sif_over_par-sif_over_par_sd)/4600,ymax=(sif_over_par+sif_over_par_sd)/4600)) +
  #scale_color_manual(values = c("#BB5566","#DDAA33","#004488","#888888")) +
  #scale_color_manual(values = c('#999999', '#FB9A29', '#EC7014', '#993404') ) +
  labs(x=bquote("Total Chlorophylls"~~~"["*mu*mol~m^-2*"]"), y=bquote("Total Xanthophylls"~~~"["*mu*mol~m^-2*"]"))+
  theme(legend.position = "none")

ggplot(merged_df,aes(x=par,y=sif_fld,color=long_desc)) +
  geom_point() +
  scale_color_manual(values=c('#999999','#E7D4E8','#C2A5CF', '#9970AB',
                              '#FEE391', '#FEC44F', '#81C4E7')) +
  theme_tsj() +
  labs(x=bquote("PAR"~~~"["~mu*mol~m^-2*s^-1~"]"), y=bquote("SIF"~~~"["~mW~m^-2*sr^-1*nm^-1~"]")) +
  geom_errorbar(aes(ymin=(sif_fld-sif_fld_sd),ymax=(sif_fld+sif_fld_sd))) +
  geom_errorbarh(aes(xmin=(par - par_sd),xmax=(par + par_sd))) +
  geom_abline( slope=.000436, intercept = 0.0625,size=.4,linetype=2) +
  geom_abline( slope=.000429, intercept = 0.0151,size=.4) +
  theme( legend.position = "none")

sub_df <- merged_df[ merged_df$long_desc == "A", ]
fit <- lm( sif_fld~par, sub_df)
summary(fit)

fit <- lm( sif_fld~par, merged_df)
summary(fit)

licor_file <- "licor_data.csv"
licor_df   <- load_licor_data( licor_file )
licor_df$tstamp <- as.POSIXct( licor_df$tstamp + 4*60*60, origin="1970-01-01", tz="UTC")
licor_df$target_type <- factor( licor_df$target_type,levels=c("Control","ABA","Drought"))

fit   <- lm(PhiPS2~A,licor_df)
coefs <- coefficients(fit)

licor_df$long_desc <- "Other"
licor_df$long_desc[ licor_df$target_type == "Drought" & licor_df$date == "20210716"] <- "Drought: Rewatered"
licor_df$long_desc[ licor_df$target_type == "Drought" & licor_df$date == "20210715"] <- "Drought: Day 2"
licor_df$long_desc[ licor_df$target_type == "ABA" & licor_df$date == "20210716"] <- "ABA: Day 3"

g <- ggplot(all_df, aes(y=PhiPS2,x=A,color=long_desc)) +
  geom_point() +
  geom_errorbarh(aes(xmin=A-A_sd,xmax=A+A_sd))+
  geom_errorbar(aes(ymin=PhiPS2-PhiPS2_sd,ymax=PhiPS2+PhiPS2_sd))+
  ylim(0,.2) +
  #scale_color_manual(values=c('#9970AB','#FEC44F','#999999')) +
  scale_color_manual(values=c('#999999','#E7D4E8','#C2A5CF', '#9970AB',
                              '#FEE391', '#FEC44F', '#81C4E7')) +
  geom_abline( intercept = 0.03917, slope = .009548,size=.4) +
  geom_abline( intercept = 0.043397, slope = .009139,size=.4,linetype=2) +
  labs(y=bquote(phi[PSII]), x=bquote(A[net]~~~'['*mu*mol~m^-2*s^-1*']'), color=NA) +
  #theme(legend.position = c(.75,.25)) +
  #theme(legend.position = "none") +
  #theme(legend.key = element_rect(fill = "white")) +
  theme(legend.background = element_rect(fill = "white")) +
  theme(legend.text = element_text(  size = rel(.8))) +
  theme(legend.title=element_blank() ) +
  theme(axis.text = element_text(  size = rel(.9))) +
  theme(axis.title = element_text(  size = rel(1.1)))
  #annotate("text",x=2.25,y=.185,label="italic(R) ^ 2 == 0.89", parse=T,size=5) 
g
ggsave("plots/A_vs_phiPSII.jpg",g,dpi=300,width=5,height=4,units="in")

foo <- all_df[ all_df$long_desc !="C4", ]
g <- ggplot(foo, aes(y=SIF_yield,x=FvFm,color=long_desc)) +
  geom_point() +
  #geom_errorbarh(aes(xmin=A-A_sd,xmax=A+A_sd))+
  #geom_errorbar(aes(ymin=PhiPS2-PhiPS2_sd,ymax=PhiPS2+PhiPS2_sd))+
  #ylim(0,.2) +
  #scale_color_manual(values=c('#9970AB','#FEC44F','#999999')) +
  scale_color_manual(values=c('#999999','#E7D4E8','#C2A5CF', '#9970AB',
                              '#FEE391', '#FEC44F', '#81C4E7')) +
  #geom_abline( intercept = 0.03917, slope = .009548,size=.4) +
  #geom_abline( intercept = 0.043397, slope = .009139,size=.4,linetype=2) +
  labs(y=bquote(SIF[yield]), x=bquote(Fv*"'"*"/"*Fm*"'"), color=NA) +
  #theme(legend.position = c(.75,.25)) +
  #theme(legend.position = "none") +
  #theme(legend.key = element_rect(fill = "white")) +
  theme(legend.background = element_rect(fill = "white")) +
  theme(legend.text = element_text(  size = rel(.8))) +
  theme(legend.title=element_blank() ) +
  theme(axis.text = element_text(  size = rel(.9))) +
  theme(axis.title = element_text(  size = rel(1.1)))
#annotate("text",x=2.25,y=.185,label="italic(R) ^ 2 == 0.89", parse=T,size=5) 

g

ggsave("plots/A_vs_phiPSII.jpg",g,dpi=300,width=5,height=4,units="in")


fit <- lm(PhiPS2~A,all_df)
coefs <- coefficients(fit)

sub_df <- all_df[ all_df$long_desc == "A", ]
fit <- lm(PhiPS2~A,sub_df)
coefs <- coefficients(fit)

fit <- lm(PhiPS2~A,licor_df, weights = 1/(coefs[2]*A_sd + PhiPS2_sd) )
summary(fit)

fit <- lm(SIF_yield~FvFm,foo)
coefs <- coefficients(fit)
fit <- lm(SIF_yield~FvFm,foo, weights = 1/(coefs[2]*FvFm_sd + SIF_yield_sd) )
summary(fit)


# ABA DOSAGE EXPERIMENT-------------------------------------------------------------------------------

dose_pal  <- c("#888888", "#58508d", "#bc5090", "#ff6361", "#ffa600")

in_file <- "aba_dose_test.csv"
df <- read.csv(in_file,stringsAsFactors = F)
df$ANET_var <- df$ANET_SD^2
df$ANET_w_var <- (df$ANET_N - 1)*df$ANET_var
df$ANET_nx <- df$ANET*df$ANET_N
df$neg <- df$ANET
df$neg[ df$Rep == "B"] <- -1*df$neg[ df$Rep == "B"]
sum_df  <- aggregate( cbind(ANET,ANET_nx,ANET_N,ANET_w_var,neg)~Day+Dose,df,"sum")
min_df  <- aggregate( ANET~Day+Dose,df,"min")
prod_df <- aggregate( ANET_N~Day+Dose,df,"prod")
gcc_df  <- aggregate( Gcc~Day+Dose,df,"mean")
sd_df   <- aggregate( Gcc~Day+Dose,df,"sd")

mdf <- sum_df
mdf$ANET <- sum_df$ANET_nx/sum_df$ANET_N
mdf$ANET_var <- sum_df$ANET_w_var/(sum_df$ANET_N - 1) #+ prod_df$ANET_N*((sum_df$neg)^2)/(sum_df$ANET_N*(sum_df$ANET_N-1))
mdf$Dose <- as.character(mdf$Dose)
mdf$ANET_sd <- sqrt(mdf$ANET_var)
gcc_df$Gcc_sd <- sd_df$Gcc
gcc_df$Dose <- as.character(gcc_df$Dose)

g <- ggplot(mdf,aes(x=Day,y=ANET,color=Dose)) +
  geom_point() +
  geom_line() +
  geom_errorbar( aes(ymin=ANET-ANET_sd,ymax=ANET+ANET_sd),width=.2,alpha=.5 ) +
  scale_color_manual(values=dose_pal) +
  scale_x_continuous(breaks=0:10) +
  scale_y_continuous(breaks=seq(-2,10,2)) +
  theme_bw() +
  labs(x="Days since Application",y=bquote(A[net]~~~'['*mu*mol~m^-2*s^-1*']'),color="Relative ABA Dosage")+
  theme(legend.position = "bottom")
ggsave("plots/fig4_a.jpg",g,dpi=300,width=6,height=3.75,units="in",type='cairo')

g <- ggplot(gcc_df,aes(x=Day,y=Gcc,color=Dose)) +
  geom_point() +
  geom_line() +
  geom_errorbar( aes(ymin=Gcc-Gcc_sd,ymax=Gcc+Gcc_sd),width=.2,alpha=.5 ) +
  scale_x_continuous(breaks=0:10) +
  scale_color_manual(values=dose_pal) +
  theme_bw() +
  labs(x="Days since Application",y="Gcc",color="Relative ABA Dosage")+
  theme(legend.position = "bottom")
ggsave("plots/fig4_b.jpg",g,dpi=300,width=6,height=3.75,units="in",type='cairo')

















